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Abstract - (Bi)multi-partite interaction patterns are commonly observed in real world systems 
which have inhibitory and excitatory couplings. We hypothesize these structural interaction pat¬ 
tern to be stable and naturally arising in the course of evolution. We demonstrate that a random 
structure evolves to the (bi)multi-partite structure by imposing stability criterion through mini¬ 
mization of the largest eigenvalue in the genetic algorithm devised on the interacting units having 
inhibitory and excitatory couplings. The evolved interaction patterns are robust against changes 
in the initial network architecture as well as fluctuations in the interaction weights. 
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c/3 Introduction: Finding mechanisms governing the evo- 
^lution of patterns in real world systems remains challeng- 
i ng in the evolutionary science [1]. Structural pattern of 
1 —la network not only influences the functional response [5], 
but also is motivated from a specific function of the cor- 
responding system [3]. Exploring evolutionary origins of 
— existing structures is important to develop an understand- 
00 ing pertaining to the hidden roles of these structures. For 
r ^ example, nestedness in ecological networks structure has 
been shown as an evolved property in co-operative inter- 
actions [3] and also in maximizing the structural stability 
of mutualistic systems EE]- 


in 


This Letter focuses on understanding the evolutionary 
origin of bipartiteness and multipartiteness in real world 
networks. A network is called bipartite (multipartite) if 
. its nodes can be divided into two (more than two) groups 
in such a manner that nodes in one group are connected 
5^ to the nodes in the other groups with no or sparse con- 




iiection existing within the same group. There are many 
natural systems which posses the bipartite (or multipar¬ 
tite) network structures. Few examples of the real world 
networks, which are build up to the current state as a re¬ 
sult of evolution and possess the bipartite (or multipartite) 
structure, are gene regulatory interaction networks con¬ 
sisting of two groups of regulating and regulated genes [7] , 
ecological food webs having different trophic levels Eli, 
El Verde rainforest m, etc. However, understanding be¬ 
hind the organization of different trophic levels is still not 
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adequate m- 

Additionally, several systems are known to be character¬ 
ized by their constituents defined by different behaviors. 
For example, according to Dale’s principle, in a neural 
system all outgoing synapses of a neuron are either in¬ 
hibitory (I) or excitatory (E) [T3] , thereby leading to vari¬ 
ous structural properties due to FE, FI and E-E couplings 
[T3l|T4]. Whereas, in ecosystem, several species show bias 
for the mixed behavior. Eor example, Apex predators have 
zero or less number of predators for their own community 
|15] . These different behaviors of nodes or interactions are 
crucial for a desired functioning or the robust functional 
performance of the underlying system [inillZ] ■ The cele¬ 
brated work by Robert May demonstrates that the largest 
real part of eigenvalues (i?max) of the corresponding Ja¬ 
cobian matrices contain information about the stability of 
the underlying systems [18) . The emergence of complex 
structural properties has been studied using Rmax which 
is considered as a measure defining fitness of the under¬ 
lying system during evolution [19] . Further, the Genetic 
algorithm (GA) is a widely used technique in optimiza¬ 
tion problems as well as in providing evolutionary mod¬ 
els for systems in many disciplines EHiin]- It has been 
shown that a stability maximization-based genetic algo¬ 
rithm leads to an emergence of the hierarchical modular¬ 
ity in a network [32]. Recently, evolution of clustering 
is demonstrated by maximizing the stability of the un¬ 
derlying network using GA m- This paper provides an 
explanation to the emergence of bipartite (multipartite) 
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Fig. 3: Schematic diagram for (a) initial population of ER 
networks and (b) finally evolved bipartite structnre (dashed 
circles). Dots and arrows represent I and E links, respectively. 


Fig. 1: (Color online) Mechanism to generate a child matrix 
from a pair of randomly selected fitter networks. The blocks 
of the parent and child matrices at a particular position are 
represented with the same color. The dimension of each square 
block is 10. Note that the GA works fine for all the block sizes 
signihcantly less than N. For larger block size, GA becomes 
meaningless as the child matrices lose variability and turn out 
to be very similar to their parent matrices. 


probability pin, type-II nodes (i.e. nodes with only I con¬ 
nections) are randomly selected out of the total N nodes 
(Eq. [T]). Rest of the nodes form a set of the type-I nodes 
(i.e. nodes with only E connections) of the behavior ma¬ 
trices. This arrangement leads to the following matrix, 
termed as behavior matrix, which has a balance of the 
inhibition and excitation |241I25) in the whole network, 


structure using GA-based minimization of Rmax- 

Theoretical Framework: We start with the P in¬ 
dependently drawn adjacency matrices, [oij] (each having 
dimension N) corresponding to Erdos-Renyi (ER) random 
networks representing the initial population. Elements in 
the corresponding adjacency matrices [oy] take value 1 
and 0 depending upon whether there exists a connection 
between ith and jth nodes which is decided with a proba¬ 
bility q. Next, we generate matrices in which each node is 



Fig. 2: (Golor online) (a) and (c) evolution of Rmax and 
fintra//inter respectively for Pin = 0.50 and connection weight 
from Eq. ©■ Minimum, average and maximum values for the 
network population are indicated by □ (red), o (green) and o 
(black) respectively, (b) and (d) depict evolution of Rmax and 
fintra //inter , respectively for connection weight with fluctua¬ 
tions as imposed by Eq. ©. Initial ER networks have (k) — 6 
with N = 100. For all the cases network population is 500. 

either excitatory (type-I) or inhibitory (type-II) through¬ 
out the evolution. This is motivated from the Jacobian 
(synaptic) matrices considered in the model of [24] , which 
incorporates Dale’s principle. Such type of matrices have 
following properties: (i) The columns are either Gaussian 
distributed numbers with negative mean or positive mean 
with certain variance, (ii) Summation of the numbers in 
each row is equal to zero, which takes into account the 
balanced condition. Such matrices lead to the underly¬ 
ing networks that are globally connected. Therefore, with 


j 1 - l/PinVj, if i G type-II 
I 1 Vj, if i G type-I 

More specifically in Eq.[TJ rows of matrices [fey] have either 
‘-I-1’ or T — 1/pin’ entries. The entries of the 6ij matrices 
are fixed in the course of evolution, demonstrating the in¬ 
variant nature of nodes and consequently the links [26) . 
We define another matrix [cy] (Eq. [2|) for assessing the 
fitness of a network, which is constructed using the ad¬ 
jacency matrices of the network and the behavior matrix 
(Eq. [IJ in the following manner: 


Cii 


6ij if Oij ^ 0 

0 otherwise 


( 2 ) 


A network corresponding to lower i?max value of an asso¬ 
ciated [cij] matrix is referred here as a fitter network. For 
the next generation, the P fitter networks of size N x N 
each are created as follows. First of all we select P/2 num¬ 
ber of distinct networks from these networks which have 
Rma,x values lesser than those of the rest of the networks in 
the present population. These selected networks are fitter 
networks and act as the parent networks for creating rest 
P/2 child networks for the next generation (Fig. [1}. For 
example, in order to generate one child (say C), a pair of 
networks (say A and B) from P/2 population are selected 
randomly. Then each block of the adjacency matrix of the 
child at a specified dimension and location (in terms of 
rows and columns), is filled with the block having the same 
position in the parent matrices, with an equal probability 
(Fig. [1}. The child matrices generated are asymmetric in 
nature. We thus take the upper triangular part of the 
child matrices and construct their adjacency matrix such 
that they are symmetric. In order to avoid trivial changes 
in the i?max values due to the change in the total number 
of connections in the network and to get a signature of 
the change in the network pattern on the i?max values, we 
preserve the average degree ((fc)) of the child matrices by 
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randomly removing or inserting connections with a small 
probability. This probability is decided by the fluctuation 
of the average degree in the crossed child network from 
the degree of the initial networks population used in the 
GA. Thus, the total number of edges remain almost same 
for the child networks. 

In order to emphasize that here changes in the value of 
Rmax are not due to the change in the average degree m 
but arising due to the changes in the interaction patterns, 
we present the results for the preserved average degree of 
the child network. Since real networks are known to be 
sparse in nature, the results are presented for (fc) = 6. 
A lesser value of (fc) leads to disconnected components, 
in turn affecting Rmax- With an increase in (k), R^^k 
increases in general. The larger (k) leads to a less varia¬ 
tion among the structural conformations of initial network 
population, which results in a lesser variation in the child 
population as well. The small amount of variation present 
in the set of initial networks population vanishes after cer¬ 
tain generations without producing sufficient structural 
changes with respect to the initial population. The evo¬ 
lution stops when there is no structural variation in the 
network population. Also bipartite structure does not ex¬ 
ist for a network having (k) > N/2. Thus, the GA ap¬ 
proach which preserves the average degree, does not lead 
to evolution of bipartite structure for very dense networks. 
Further, with larger N, the values of i?max may vary but 
the overall behavior that Rmax minimization leads to evo¬ 
lution of the bipartiteness is expected to remain same, 
provided the initial networks consists of different types 
of nodes. The problem encountered on considering larger 
network size [22] , would be requirement of memory alloca¬ 
tion for a population of P networks in GA amounting to 
size of the 1-d array, which is computationally not 
favorable for large N. 

Further, the total number of connections within the 
same population (within type-I or within type-II) is de¬ 
noted by /intra and the total number of connections among 
the different populations (between type-I and type-II) is 
denoted by /inter- The ratio of /intra and /inter provides 
a measure for the bipartivity [2811^ . The zero value of 



Fig. 4: (Color online) Amax during the evolution of matrices 
[cij] for children and parents population. Square (black) and 
circle (red) represent maxima and minima of Amax values of the 
[cij] matrices for the parents population respectively and dia¬ 
mond (green) and star (blue) are maxima and minima of Amax 
values of [cq] matrices for the children population respectively. 


fintra/ finter iudicatcs a bipartite structure. However the 
converse is not true, a bipartite structure might also arise 
when in a network, there are two groups of nodes such that 
one group (having both type-I and type-II category nodes) 
forms all the connections with the other group (again hav¬ 
ing both type-I and type-II category nodes) and no con¬ 
nection within the group. Gonnections between type-I 
and type-II nodes of individual groups would contribute to 
f intra leading to non-zero f intra /f inter value. Note that, 
a bipartite structure may exist in a system having only 
single behavior. For example, a system in which plants 
and pollinators form two different communities exhibiting 
the bipartite structure with only type-I couplings pgUTT]. 
This situation does not arise in our analysis as we consider 
different classes of nodes and bipartiteness is confirmed 
only when the value of f intra /f inter is zero. 

Emergence of bipartite structure: Starting with 
sparse ER random networks as an initial population, as 
evolution progresses through GA (Fig. [2|), the minimiza¬ 
tion of Amax occurs. In the initial population, the mean of 
fintra/finter is close to One, however maxima and minima 
are different as expected from the various realizations of 
the ER networks. As evolution progresses, the mean and 
minimum values of fintra/finter converge towards zero, 
however maxima still exhibits a small separation from the 
zero (Fig. [2]). This indicates that most of the networks 
in the evolved population attain the bipartite structure in 
which no or a very few connections exist within type-I or 
within type-II category nodes. Since, the GA minimizes 
Amax values of the population, the convergence of the max¬ 
ima towards zero is stopped when the minima closely co¬ 
incides with the mean value. What follows that the initial 
population consisting of random architecture (Fig. EJa)) 
converges to the bipartite structure (Fig. EKb)) through 
the stability maximization. 

The emergence of the bipartite structure directly follows 
from the spectral properties of antisymmetric matrices en¬ 
tailing all imaginary eigenvalues. The [cij] matrix, for a 
ideal bipartite structure without any connection within 
type-I or within type-II category nodes, corresponds to a 
antisymmetric matrix, i.e. [cy = - cji]. The Amax of [cij] for 
the initially considered ER networks, generated using con¬ 
nection probability q, takes positive values given by y/Nq 
for Pin = 0.50 [20]. Additionally, the correlation in the el¬ 
ements of c matrix, given as r = Y/^,j=i CijCji/ GjGj 

for the matrices defined by Eq. (T] takes value 

^ — (/intra /inter)/(/intra “b /inter) (3) 

for p„ = 0.5. This value decreases with fintra/finter, 
further leading to a decrease in Amax because the spectral 
distribution of [cji] matrices are of elliptical shape with 
its axis lying on the real line decreasing around a fixed 
center [31]. Since we minimize Amax, the GA successfully 
and smoothly evolves a structure in which fintra/finter 
has very small or zero value. 

The successful performance of the GA can be explained 
by considering two different classes of network population. 


P-3 




Sarika Jalan^ and Sanjiv K. Dwivedi 


Let the children and their parent populations form two 
different classes of networks in each generation and the 
-Rmax of [cij] associated with both the classes, decrease 
during the evolution. The maxima and minima of i?max 
for the parent [cy] decrease as per the GA minimization 
(Fig.|3]). The maxima and minima of i?max for children [cij] 
also decrease with the evolution. It happens so that the 
maxima of i?max values for children population is greater 
than the maxima of those for the parents population. The 
hamming distance between the adjacency matrices of net¬ 
works provide a method to measure this dissimilarity |32j . 
What follows that, whenever parents exhibit dissimilarity, 
i.e. Hamming distance is large, their crossed child is also 
likely to have a high deviation in its structural properties 
which in turn leads to a large deviation of i?max of the 
child from its parent. 

Furthermore, the minima value for children and parents 
populations may lead to the following two cases. In the 
first case, the minima value for children is higher than 
the minima value for the parent population. Hence, the 
minima value for parent population in the subsequent gen¬ 
erations does not decrease until the minima values for the 
children become lower than the minima value of their par¬ 
ents population. Whereas in the second case, the minima 
values for parents are higher than those of child, hence 
parents, which produce child for the next generation, are 
selected from the child in the current generation leading to 
a decrease in the minima of the parents for the next gener¬ 
ation. This leads to the value of i?max decreasing in each 
generation which indicates that as evolution progresses, 
the evolved feature (i.e. a lesser value of fintra/finter) 
is transferred to the child. It occurs due to the proce¬ 
dure of block selection for child matrices as represented in 
the schematic diagram (Fig. [ij. The GA-based optimiza¬ 
tion method converts the entries of the adjacency matri¬ 
ces corresponding to the connections within type-I nodes 
and the connections within type-H nodes to E-I couplings, 
thus leading to a decrease in fintra//inter value. The de¬ 
creased values of fintra/finter of parent matrices give rise 
to a chance to produce a child with a further decreased 
value. And whenever parents produce such type of child 
network, they possess lower iJmax values which are selected 
for the next generation. In other words, fitter children are 
produced by the fitter parents. 

Robustness against random fluctuations: To 
demonstrate the robustness of the results against changes 
in coupling strength, we introduce random fluctuations 
[I7l[33], which can be imposed in the fitness matrix as 
follows 

fx6ij ifaijyfO 

1 0 otherwise Oij =0. 

where X is a uniform random variable lying between 0 and 
1. As illustrated in Fig. [21 the bipartite structure evolved 
through evolution remains unaffected by the random fluc¬ 
tuation in coupling strength. The only implication of hav¬ 
ing the random fluctuations is that the evolutionary pro- 



Fig. 5: (Color online) (a) Average Amax and (b) r = 
{fintra//inter) for both the initial (circles), and the evolved 
(stars) networks as a function of d. The system is evolved for 
5000 iterations consisting of 500 networks population for GA. 
The average values, for the evolved networks, is calculated for 
all the population of last 1000 generations. Initial ER networks 
have (k) = 6 and N = 100. 


cess becomes bit slower. The minima, maxima and mean 
of i?max decrease during the course of evolution with the 
mean being bit closer to the minima rather than the max¬ 
ima (Fig. [2]). After few generations, three values become 
saturated i.e. parallel to the generation axis. The minima 
does not converge to zero which is unlikely for no fluctu¬ 
ations case. The structural changes during the evolution 
lead to the decrease in the minima, mean and maxima, 
whereas fintra/finter still remains close to zero similar to 
the no fluctuation case. 

Further, in order to analyze the effect of fluctuations in 
the coupling nature of the predefined coupling behaviour 
in matrix {\hij\)^ we take the range of uniform random 
variable in the negative side as well using the following 
equation; 

Y = X-d (5) 

where d is a number lying between 0 and 0.5. Therefore, 
fitness matrix [cij] would be defined as; 

Ci-= (6) 

|0 otherwise ay = 0. 

We run the simulations for various d values in which fit¬ 
ness matrix is defined by Eq. [51 For d = 0.5, matrices 
completely destroy the impact of predefined information of 
the behavior matrix([6ij]) during the evolution from Eq.|5] 
What follows that we do not get any structural changes 
during the evolution and hence there is no decrease in 
Rmax- However, for d < 0.5, the impact of behavior ma¬ 
trix on [cij] matrix is dependent on the value of d. Fig.jSKa) 
depicts that the efficiency of optimization depicted by dif¬ 
ference in the value of Rmax of the evolved and the initial 
networks, increases with a decrease in the d values. For 
d = 0.5, the rate of selection and rejection of the networks 
in a given generation becomes equal and very fast due 
to equal strength of the fluctuations in the positive and 
negative side of the entries of the fitness matrix [c^] for 
any predefined coupling behaviour {[bij]). Consequently, 
GA fails to minimize Rmax value and as a result no vari¬ 
ation in values of fintra/finter E observed at this value of 
d (Fig. |Si;a)-(b)). As d decreases, fintra/finter decreases 


p-4 







Emergence of (bi)multi-partiteness in networks having inhibitory and excitatory couplings 



Fig. 6: (Color online) Evolution of bipartite structure for sys¬ 
tem having both types of the nodes, same as that of Fig. [2] 
for the initial SF networks with {k) = 6 and N = 100. The 
network population used in GA is 500. (a) and (c) present the 
evolution of i?max and fintra /finter, respectively, for fitness 
matrix given as Eq. 2, whereas (b) and (d) present the evolu¬ 
tion of these quantities for fitness matrix being represented by 
Eq. 4. 


with a faster rate in the evolved networks (Fig. [SKa)). Af¬ 
ter a certain value of d the rate of decrement becomes 
slower. Therefore, there exists a large regime of d in which 
the evolved networks exhibit the bipartite or close to the 
bipartite structure (Fig. [SKb)). 

In order to demonstrate the robustness of results against 
changes in the initial population networks, we study the 
evolution of scale-free networks. These networks are gen¬ 
erated using Albert-Barabasi model [31] for both the cases 
of with and without fluctuations in the coupling strength. 
Since the scale-free networks have randomness arising due 
to the algorithm, it leads to the appearance of structural 
variations in different realizations. The given position 
of blocks in associated adjacency matrices have sufficient 
variations, as a result the crossed child matrices in GA at¬ 
tain conformation structure leading to bipartite arrange¬ 
ment. For coupling with random fluctuations, there is 
again an emergence of the bipartite structure (Fig. |6|). 

We have also analyzed the case when the directed net¬ 
works form the initial population for the GA. For random 
directed networks, there would be a very less probabil¬ 
ity of finding un-directional {Aij=Aji) coupling and this 
probability is dependent on the size and average degree 
of the network. However, in real world networks, such 
un-directional couplings are found with very high proba¬ 
bility [5S] . For initial networks being directed, the evolved 
networks exhibit bipartiteness as exhibited by undirected 
initial networks as long as the child networks of each gen¬ 
eration preserve directionally possessed by the parent net¬ 
works. The evolution to the bipartiteness as a result of sta¬ 
bility maximization directly follows from the discussions 
around Eq. 3. The only difference is that it requires more 
number of generations to obtain the evolved networks as 
compared to the undirected case. 

Addition of type III nodes and emergence of 
multipartite trophic networks: The results presented 
above restrict nodes with only inhibitory and excitatory 



Fig. 7: (Color online) (a) Average r for the evolved bipartite 
networks as a function of p. The population at each step con¬ 
sists of 500 networks, (b) Change in the average value of total 
number of edges among type-I and type-II nodes over the net¬ 
work population depicted by Qintra with the evolution. Zero 
value of Qintra indicates formation of different trophic levels 
consisting type-I and type-II nodes each, (c) Average behavior 
of i?max over network population with the evolution. Calcula¬ 
tions are done for the size of network being 100 and average 
degree 6 with type-I, type-II and type-III being 20, 20 and 
60, respectively. Network population used for (b)-(c) is 1000. 
For all the cases, [ci^j matrices are defined by Eq. 4. (d) Ini¬ 
tial population of ER networks when all three types of nodes 
(type-I, type-II and type-III) coexist (e) finally evolved struc¬ 
ture consisting different trophic levels (dashed circles). Dots 
and arrows represent I and E links, respectively. 


connections, which is motivated from neural models |24] . 
If we relax this criterion, nodes may have inhibitory as 
well as excitatory connections leading to characterization 
of nodes as type-III (Fig. IZKd)). We investigate the evo¬ 
lution of such type of networks. Starting with a matrix 
consisting of 1 and —1 entries randomly distributed, which 
corresponds to all nodes being type-III, there may be some 
rows having more 1 and some rows having more —1 en¬ 
tries. Now we select a row having more 1 entries and 
with probability p convert —1 entries into 1 entries. Simi¬ 
larly, for rows having more —1 entries, with probability p 
we convert 1 entries into —1 entries. The value of p be¬ 
ing zero corresponds to the original matrix which we have 
started with denoting all type-III category nodes, whereas 
p = 1 leads to the situation where those rows having more 
1 entries in the initial matrix, have all their —1 entries 
converted into 1 corresponding to type-I nodes and those 
rows having more — 1 entries in the initial matrix, have all 
their 1 entries converted into — 1 corresponding to type-II 
nodes. We find that with increase in p, there is increase 
in bipartiteness in the evolved networks (Fig. EKa)). 

In the above case, for all values of p either we have only 
type-III nodes or type-I and type-II nodes with a feeble 
probability of co-existence of all three types. Whereas 
in a realistic situation, for instance in ecological systems, 
there exist few species which are only prey (type-I), few 
others which are only predator (type-II) and some more 
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species which are prey to some species and act as predators 
for others (type-III), for example, the herbivores in the 
grassland food web m- In order to mimic this realistic 
situation which is a mixture of all three types, we consider 
some of the nodes as type-I, some of the nodes as type-II 
and rest as type-III. Again type-I nodes will lead to all 1 
entries and type-II nodes will lead to all — I entries in the 
behavior matrix. In rows corresponding to type-III nodes, 
1 and — I entries are introduced with the equal probability. 

We find that with an increase in the generations, connec¬ 
tions within type-I nodes and connections within type-II 
nodes decrease and converge towards zero (Fig. [T^b) and 
Fig. EKe)). This leads to a situation when type-I nodes 
have type-II or type-III neighbors and type-II nodes have 
type-I or type-III neighbors. GA minimizes Rmux with 
relatively faster rate in the beginning of the evolution 
and as the evolution progresses the rate becomes slower 
(Fig. EKc)) due to a steep decrease in the value of r ^5] . 
Further, number of the type-III nodes, which are neighbors 
of the type-I nodes but are not neighbors of the type-II 
or vice versa, are maintained with random fluctuations in 
the course of evolution. Additionally, with the evolution, 
the number of connections among such type-III nodes de¬ 
creases. The nodes of type-I and type-II are not connected 
within their categories and thus form two distinct groups 
(Fig. EKe)). Also, the type-III nodes, which are neigh¬ 
bors of type-I (type-II) nodes but are not the neighbors 
of type-II (type-I) nodes (Fig. EKe)), form other group(s). 
This organization resembles the various trophic levels of 
food-web where the species belonging to same trophic level 
do not have any interaction among themselves [8l[9] 

Conclusion: To conclude, we present an evolutionary 
origin of the emergence of bipartiteness in interaction pat¬ 
terns. The essence of the method lies in the property that 
individuals get classified in to different classes in which 
they have hxed invariant behaviors. The Rmax values are 
dependent on network architecture as well as fluctuations 
in coupling strength. Since Rmax quantifies the stabil¬ 
ity of a system [T5], there arises a chance of this value 
to increase, making a system unstable. Additionally, an 
increase in the system size may also lead to an increase 
in Rmax implying further that a system with larger size 
ideally should not exist [H]. However, our results sug¬ 
gest that despite larger system size and high disorder in 
coupling strength, a balance of I and E limits Rmax in the 
evolved networks, and as a result the fluctuations may not 
destabilize the system. This result provides an insight in 
to why despite possessing high degree of disorder ecosys¬ 
tems are robust, which can be attributed to the existence 
of various trophic levels of species preserving stability of 
the system. 

The method considered here may be useful in optimizing 
various other measures related with structural as well as 
spectral behavior [5511571 . The framework can be extended 
further to derive an evolutionary understanding of how a 
systems’s function gets affected by underlying structural 
patterns. For instance, assortative-disassortative mixing 


interaction patterns are seen in real world networks [3], 
evolution of which can be addressed from the stability 
point of view. 
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